3次元点集合データ(ポイントクラウド)の解析

この文章で解説するのは

  1. 3次元の点データからパーシステント図を計算する
  2. その図をパラメータを変えながら可視化する
  3. 各birth/death timeを抽出する
  4. 逆解析を行う

の4つです。ここまではおよそどのようなポイントクラウドデータでも 共通ですので、ここまでできるようになると良いでしょう。

パーシステント図の計算

対称となるデータは pointcloud.txt というファイルです。これを解析してみましょう。

まず最初に、必要なライブラリをインポートします。

In [1]:
import homcloud.interface as hc  # HomCloudのインターフェス
import homcloud.plotly_3d as p3d  # 3次元可視化用
import plotly.graph_objects as go  # これも3次元可視化用
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

データを np.loadtxt で読みこみます。

In [2]:
pointcloud = np.loadtxt("pointcloud.txt")

まずは pointcloud の情報を調べてみましょう．

In [3]:
pointcloud.shape
Out[3]:
(1000, 3)
In [4]:
np.min(pointcloud, axis=0), np.max(pointcloud, axis=0)
Out[4]:
(array([0.00026274, 0.00100671, 0.00020878]),
 array([0.99997803, 0.9993682 , 0.99825316]))

データは 1000x3 の配列です．これは3次元の点が1000個あることを意味します． またX，Y，Z座標の最大最小を見るとデータは$[0, 1]\times[0, 1]\times[0, 1]$という立方体に分布していることがわかります． この点を可視化してみます．

In [5]:
go.Figure(
    data=[p3d.PointCloud(pointcloud, color="black")],
    layout=dict(scene=dict(xaxis=dict(visible=False), yaxis=dict(visible=False), zaxis=dict(visible=False)))
)

可視化には plotly というライブラリを使っており，ブラウザ上のjupyter notebookで直接表示できます． マウスの左ボタンドラッグで回転，右ボタンドラッグで平行移動，ホイール回転で拡大縮小ができます．

dict(xaxis=dict(visible=False), yaxis=dict(visible=False), zaxis=dict(visible=False)) というのを略記するために p3d.SimpleScene という関数が使えます．

In [6]:
# 上と同じ
go.Figure(
    data=[p3d.PointCloud(pointcloud, color="black")],
    layout=dict(scene=p3d.SimpleScene())
)

データからパーシステント図を計算します。hc.PDList.from_alpha_filtrationという静的メソッドを使います。 save_boundary_map は後で optimal volume というものを計算するのに必要です。

In [7]:
hc.PDList.from_alpha_filtration(pointcloud, 
                                save_to="pointcloud.pdgm",
                                save_boundary_map=True)
Out[7]:
PDList(path=pointcloud.pdgm)

すると pointcloud.pdgm というファイルが生成されます。これが パーシステント図の情報を収めたファイルです。

注意

HomCloud 3.0 以降では .pdgm という拡張子を推奨しています．

パーシステント図の可視化

次に、ここからパーシステント図をプロットしてみましょう。

まず、ファイルを読み込みます。

In [8]:
pdlist = hc.PDList("pointcloud.pdgm")

このオブジェクトはすべての次数のパーシステント図を保持している(0次から2次まで)ので、1次のパーシステント図だけ取り出します。

すべての次数を保持しているので、このクラス自体の名前はPDListという名前です。

In [9]:
pd1 = pdlist.dth_diagram(1)

このメソッドの返り値は PD class のインスタンスです。

ヒストグラムを構築し、パーシステント図をプロットします。 histogramメソッドでヒストグラムを構築し、plotでそれをプロットします。

In [10]:
pd1.histogram().plot()

左下のほうに何かちょっと見えてあとはかなりぼんやりしています。というのは このデータは左下のほうにデータ(birth-death pair)が偏っているからです。 そこで色付けを log scale にしましょう。

In [11]:
pd1.histogram().plot(colorbar={"type": "log"})

基本的にはパーシステント図の点は対角線から離れるほど「意味のある」リング構造と対応し、 またY軸が大きい値になるほど大きなリング構造を表しています。 つまり(0.0025, 0.008)のあたりにある点がリング構造のなかで最もちゃんとしたリング構造を持っているものと対応していると言えそうです。

ここで注意しておくと、パーシステント図のX、Y軸は点に貼り付ける球の半径と対応していると よくある教科書には書かれています。HomCloudでは実はこれは半径の2乗が使われています。 つまり$\sqrt{0.0025}=0.5$と$\sqrt{0.008} \approx 0.09$ が実際の半径になります。これは主には内部のアルゴリズムの都合に よるものですが、各点に重みづけをしたときにはこのほうが自然に見えるという事情もあり、 HomCloudでは2乗の値が使われます。これを止めたいときはパーシステント図をhc.PDList.from_alpha_filtration で計算するときに no_squared=True という引数を付けると 半径パラメータそのものがX、Y軸に現れます。

さて、次にbirth-death pairが集中している 左下の 0.0〜0.01 のあたりを拡大して調べてみましょう。デフォルトでは ヒストグラムを構築するときにすべてのbirth-death pairが含まれるような 範囲を自動的に切り取ります。これを変更するには histogram メソッドのx_range引数（第1引数）を指定します。

In [12]:
pd1.histogram((0, 0.01)).plot(colorbar={"type": "log"})

ヒストグラムの解像度を変更するには histogram メソッドの x_bins引数（第2引数）を指定します。デフォルトでは 128x128でヒストグラムを描きますが、もっと細かく256x256にしてみましょう。

In [13]:
pd1.histogram((0, 0.01), 256).plot(colorbar={"type": "log"})

これらの図は matplotlib の機能を使って描画されているため、 matplotlib.pyplot.savefigで保存することができます。以下のようにすると、pointcloud-pd1.pngに図が保存されます。

In [14]:
pd1.histogram((0, 0.01), 256).plot(colorbar={"type": "log"})
plt.savefig("pointcloud-pd1.png")

2次元のPDも同様に可視化しましょう．

In [15]:
pd2 = pdlist.dth_diagram(2)
In [16]:
pd2.histogram().plot(colorbar={"type": "log"})

こちらは death (Y軸) が 0.0125 より大きいあたりにぱらぱらと birth-death pair が分布しているのが特徴的です．

数値データの読み出し

pd1 や pd2 が指すオブジェクトの births、deaths属性を調べることで、birth time、death timeがわかります。

In [17]:
pd1.births
Out[17]:
array([0.00140906, 0.00139591, 0.00134116, ..., 0.01572219, 0.01709659,
       0.01871588])
In [18]:
pd1.deaths
Out[18]:
array([0.0015926 , 0.00164399, 0.00170775, ..., 0.0157662 , 0.01717465,
       0.01880101])

そこで、この2つの差を見ることでlifetimeが計算できます。

In [19]:
pd1.deaths - pd1.births
Out[19]:
array([1.83534435e-04, 2.48078852e-04, 3.66588805e-04, ...,
       4.40061409e-05, 7.80638251e-05, 8.51306753e-05])

lifetime のヒストグラムを表示してみましょう。

In [20]:
plt.hist(pd1.deaths - pd1.births, bins=100);

行末の ; としているのは plt.hist の返り値を無視するためのトリックです。大半の lifetime は非常に小さいことがわかります。実際対角線のそばに多くの birth-death pair がノイズ的な情報として分布していることはよくあります．

逆解析について

パーシステント図の個々の点は何らかの意味で元のポイントクラウドのリング構造、空隙構造と対応しているはずです。しかしそれがどのようなものであるのかを特定するのは実は そんなに簡単ではないです。このような解析を逆解析と呼びます。

ここでは HomCloud の強力な逆解析ツールである、optimal volumeを使いましょう。とりあえず (0.0025, 0.008) 付近にある birth-death pair の optimal volume を調べます。 optimal volume について詳しくは、大林の論文を参考にしてください。

pd.nearest_pair_to で (0.0025, 0.008) に一番近い birth-death pair を検索することができます。

In [21]:
pair = pd1.nearest_pair_to(0.0025, 0.008)

このオブジェクトは Pair クラスのオブジェクトで、birth timeやdeath timeといった情報を保持しています。

In [22]:
pair
Out[22]:
Pair(0.0025623095880009357, 0.008135892057712315)

以下のようにすると optimal volume が計算できます。

In [23]:
optimal_volume = pair.optimal_volume()

boundary_points メソッドで、対応するリング構造に含まれる点の座標が得られます。

In [24]:
optimal_volume.boundary_points()
Out[24]:
[[0.24041004740553085, 0.3620263927982651, 0.3938519814000714],
 [0.07850470821281752, 0.31832538225761675, 0.2726051865809148],
 [0.26866264289512143, 0.5145750814813962, 0.5502132034822121],
 [0.1255383953826733, 0.0800900311964785, 0.7190957020829821],
 [0.1252375238953729, 0.20843145021582987, 0.37381081780757186],
 [0.27896903483481006, 0.2900038108788685, 0.3853792732744701],
 [0.11968719444459341, 0.46213675229879714, 0.6086047255929676],
 [0.10439142042970573, 0.12004955606294554, 0.3999294686167837],
 [0.07478379187666995, 0.1318144315592794, 0.6167927456566142],
 [0.241996106473225, 0.4429189262369284, 0.4639170406136188],
 [0.2422121988843836, 0.4155487622456445, 0.7075748545403898],
 [0.33836973880327526, 0.5114239253242512, 0.7178588306122617],
 [0.19443369660602383, 0.2643008696498902, 0.7685550494906768],
 [0.06800441912888222, 0.1305283961003515, 0.5058547294424975],
 [0.33213749494421396, 0.4315334142137247, 0.7323967387271755],
 [0.41898602615433533, 0.44362916692289767, 0.3710002303374377],
 [0.12195279248515212, 0.003239458152972219, 0.6061209669659768],
 [0.21608561791517278, 0.2918069725141198, 0.3219258262415402],
 [0.1894939883612946, 0.03095999913576697, 0.6762565272706147],
 [0.0708147888171643, 0.05566245889432686, 0.6142169362499716],
 [0.39171929078549805, 0.4405548898038244, 0.288435385498083],
 [0.006270880292196468, 0.10979896194917471, 0.5705032405071087],
 [0.3166872139448015, 0.38348372951946463, 0.33415602643779374],
 [0.3806980910779977, 0.4089441782551778, 0.43539259723846135],
 [0.3171296008248521, 0.531306547414789, 0.6355063863352421],
 [0.15814254898557367, 0.44579878013749064, 0.4799628611505681],
 [0.17887993848973416, 0.42554668220768377, 0.6217943706156428],
 [0.09845847975583133, 0.26404248130005237, 0.7778474474038211],
 [0.3373945662185728, 0.46415106121781036, 0.4856613355037861],
 [0.08346167978435703, 0.26153343284444097, 0.3241366037151565],
 [0.03223409787023468, 0.08546385914763066, 0.4515319922985789],
 [0.08561651425017891, 0.42297356026275923, 0.528205350556188],
 [0.12886238414379725, 0.17512618876757002, 0.7438074952407614],
 [0.16660019888277078, 0.33716865816694763, 0.30111066188327684],
 [0.18595842828877163, 0.502526980520708, 0.6007699835223476],
 [0.32035651808142795, 0.4751033316678963, 0.34029734006837487],
 [0.18678772331743787, 0.34044366488617717, 0.7226556446102448]]

このリング構造は39点からなっていることがわかります。boundaryメソッドを使うとどの点とどの点がつながっているか等もわかります。

このoptimal volumeを3次元可視化してみましょう。

In [25]:
go.Figure(data=[optimal_volume.to_plotly3d()], layout=dict(scene=p3d.SimpleScene()))

リングが表示されます。

次のようにするとポイントクラウドを重ねて表示することもできます．

In [26]:
go.Figure(data=[
    optimal_volume.to_plotly3d(width=4, name="Optimal Volume"), 
    p3d.PointCloud(pointcloud, color="black", name="Pointcloud")
], layout=dict(scene=p3d.SimpleScene()))

to_plotly_meshを使うとループの内側の面を貼るような可視化も可能です．見やすい可視化を工夫してみてください． fig.update_traces(opacity=0.5, selector=dict(name="Optimal Volume")) は面の部分の不透明度を50%にしています．

In [27]:
fig = go.Figure(data=[
    optimal_volume.to_plotly3d(width=4, name="Optimal Volume outline"), 
    optimal_volume.to_plotly3d_mesh(color="red", name="Optimal Volume"), 
    p3d.PointCloud(pointcloud, color="black", name="Pointcloud")
], layout=dict(scene=p3d.SimpleScene()))
fig.update_traces(opacity=0.5, selector=dict(name="Optimal Volume"))

次のようにするとブラウザの別のタブで開きます．最後の行のfig.show(renderer="browser")が鍵です．

In [28]:
fig = go.Figure(data=[
    optimal_volume.to_plotly3d(width=4, name="Optimal Volume outline"), 
    optimal_volume.to_plotly3d_mesh(color="red", name="Optimal Volume"), 
    p3d.PointCloud(pointcloud, color="black", name="Pointcloud")
], layout=dict(scene=p3d.SimpleScene()))
fig.update_traces(opacity=0.5, selector=dict(name="Optimal Volume"))
fig.show(renderer="browser")

name=...という引数で各オブジェクトに名前を付けたり，width=4という引数で線の太さを変更したりできます．

こんどは複数の点のoptimal volumeを同時に可視化してみます．birth < 0.004, death > 0.006 くらいの領域の birth-death pair を 以下のようにして取り出してきます．

In [29]:
pairs = [pair for pair in pd1.pairs() if pair.birth_time() < 0.004 and pair.death_time() > 0.006]

取り出してきた birth-death pair は以下のようにして数えることができます．

In [30]:
len(pairs)
Out[30]:
54

では optimal volume を計算して，可視化してみましょう．54個もあるので計算には結構時間がかかります．そこで cutoff_radiusというのを指定します．これは計算ターゲットが存在する領域の大きさを指定することで計算時間を 削減するための仕組みです．だいたいこの球の内側にあるだろう，という領域の半径を指定します． この半径を予想するのは簡単ではありませんがここではデータ全体の座標が0〜1の間なので適当に1/3くらいの 0.3でいいかな，と予想してやりましょう．

In [31]:
%%time
optimal_volumes = [pair.optimal_volume(cutoff_radius=0.3) for pair in pairs]

CPU times: user 10.7 s, sys: 369 ms, total: 11 s
Wall time: 12.8 s
In [32]:
go.Figure(data=[
    v.to_plotly3d(width=4, name=f"({v.birth_time():.5f}, {v.death_time():.5f})") for v in optimal_volumes
] + [
    p3d.PointCloud(pointcloud, color="black", name="Pointcloud")
], layout=dict(scene=p3d.SimpleScene()))

2次元のPDの逆解析について

空洞のほうも同様に可視化してみましょう．birth > 0.0125 かつ death - birth > 0.0005 などがおもしろそうです．

In [33]:
pairs = [pair for pair in pd2.pairs() if pair.birth_time() > 0.0125 and pair.lifetime() > 0.0005]
In [34]:
optimal_volumes = [pair.optimal_volume() for pair in pairs]
In [35]:
go.Figure(data=[
    v.to_plotly3d(width=4, name=f"({v.birth_time():.5f}, {v.death_time():.5f})") for v in optimal_volumes
] + [
    p3d.PointCloud(pointcloud, color="black", name="Pointcloud")
], layout=dict(scene=p3d.SimpleScene()))

次のようにすると空隙部分の体積を可視化することもできます．update_tracesについては https://plotly.com/python/creating-and-updating-figures/ などを参考にしてください．

In [36]:
fig = go.Figure(data=[
    v.to_plotly3d_mesh(name=f"({v.birth_time():.5f}, {v.death_time():.5f})") for v in optimal_volumes
] + [
    p3d.PointCloud(pointcloud, color="black", name="Pointcloud")
], layout=dict(scene=p3d.SimpleScene()))
fig.update_traces(opacity=0.5, selector=dict(type="mesh3d"))